% ***************************************************************************************
% m_calibration_noaggfund_1950.m 
% ***************************************************************************************
clear all; 
close all;
format shortG;
% ==================== LOAD FUNDAMENTAL AND DATA  =========================
load parameters.mat;                     % basic parameters 
load calib_productivity_amenity.mat;     % inverted fundamentals every 5 years from 1955
load calib_setup.mat;                    % Info data every 5 years from 1945
load commuting_cost.mat;  
load calib_auxiliary.mat;                % Info used in this calibration 

Rdata45 = Rdata_all(:,1); Rdata50 = Rdata_all(:,2); 
Ldata45 = Ldata_all(:,1); Ldata50 = Ldata_all(:,2); 
N = size(Ldata_all,1); 

% =========== COMPUTE LOWER BOUNDS for Theta 1950 & SET Theta 1950 ======== 
Rthetalb = (Rdata50 - Rdata45 < 0);   Rthetalb = (1 - Rdata50./Rdata45).*Rthetalb;  Rthetalb = max(Rthetalb); 
Lthetalb = (Ldata50 - Ldata45 < 0);   Lthetalb = (1 - Ldata50./Ldata45).*Lthetalb;  Lthetalb = max(Lthetalb);
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;
theta50 = 0.90;
if min(theta50 - thetalb)<0
    disp('>>> CHECK THETA VALUES (LOWER BOUND) ! <<< ');   pause
end 
theta55 = thetavec(1); 

% ===========
l_ave_b = sum(log(bb),2)./size(bb,2); % average amenities over time (1955-75) 
l_ave_a = sum(log(aa),2)./size(aa,2); % average productivity over time (1955-75)

%  ================ INVERSION OF (X,Y) FOR 1950 =============== 
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
X0=ones(N,1);  
Y0=ones(N,1);  
Rsyst = @(x)fReval2(x,K50,Q50,Rdata45,Rdata50,Ldata45,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);
Lsyst = @(y)fLeval2(y,K50,Q50,Rdata45,Rdata50,Ldata45,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);  
[x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0);   
[y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); 
X50=abs(x_fsolve); X50=X50./geomean(X50);
Y50=abs(y_fsolve); Y50=Y50./geomean(Y50);
X50 = (X50./geomean(X50)).*exp(mean(l_ave_b)); 
Y50 = (Y50./geomean(Y50)).*exp(mean(l_ave_a)); 

% =========== SAVE RESULTS ==============
result=array2table([ID, X50, Y50, Rdata45, Rdata50, Ldata45,Ldata50, ...
   l_ave_a, l_ave_b, cbddist, Area],'VariableNames',{'geocode','X50noagg','Y50noagg',...
    'pop45','pop50','emp45','emp50','loga5575','logb5575','cbddist', 'area'}); 
writetable(result,'../../output/quant/output_noaggfund_50.csv');  

% ************************************************************************
% ************************************************************************
% **************************  FUNCTIONS **********************************
% ************************************************************************
% ************************************************************************
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) 

N = size(DRpre,1); 
X=abs(X);  
X=X./geomean(X);
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)

N = size(DRpre,1); 
Y=abs(Y); 
Y=Y./geomean(Y);
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end
